/*
 * 	This file is a modification of phylip's drawtree.c. It is also part 
 *  of MVD_GUI.
 *
 *  MVD_GUI is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  MVD_GUI is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with MVD_GUI.  If not, see <http://www.gnu.org/licenses/>.
 *  This version copyright 2010 Desmond Schmidt
 */
/* version 3.6. (c) Copyright 1993-2004 by the University of Washington.
   Written by Joseph Felsenstein, Akiko Fuseki, Sean Lamont, and Andrew Keeffe.
   Permission is granted to copy and use this program provided no fee is
   charged for it and provided that this copyright notice is not removed. */

#include <math.h>
#include <syslog.h>
#include <stdio.h>
#include <stdlib.h>
#include "phylip.h"
#include "dist.h"
#include "fitch.h"

extern FILE *infile;
extern FILE *outtree;
/**
 * Write out file with representation of final tree
 * used in fitch & neighbor
 */
void treeout(node *p, int *col, double m, boolean njoin, node *start)
{
  	int i=0, n=0, w=0;
	char c;
	double x=0.0;

	if (p->tip)
	{
		n = 0;
		for (i = 1; i <= nmlngth; i++)
		{
			if (nayme[p->index - 1][i - 1] != ' ')
				n = i;
		}
		for (i = 0; i < n; i++)
		{
			c = nayme[p->index - 1][i];
			if (c == ' ')
				c = '_';
			putc(c, outtree);
		}
		*col += n;
	}
	else
	{
		putc('(', outtree);
		(*col)++;
		treeout(p->next->back, col, m, njoin, start);
		putc(',', outtree);
		(*col)++;
		if (*col > 55) {
			putc('\n', outtree);
			*col = 0;
		}
		treeout(p->next->next->back, col, m, njoin, start);
		if (p == start && njoin)
		{
			putc(',', outtree);
			treeout(p->back, col, m, njoin, start);
		}
		putc(')', outtree);
		(*col)++;
	}
	x = p->v;
	if (x > 0.0)
		w = (long)(m * log(x));
	else if (x == 0.0)
		w = 0;
	else
		w = (long)(m * log(-x)) + 1;
	if (w < 0)
		w = 0;
	if (p == start)
		fprintf(outtree, ";\n");
	else
	{
		fprintf(outtree, ":%*.5f", (int) w + 7, x);
		*col += w + 8;
	}
}
void freetree(pointptr *treenode, int nonodes)
{
  	int i;
  	node *p, *q;

  	for (i = 0; i < spp; i++)
    	free((*treenode)[i]);
  	for (i = spp; i < nonodes; i++) 
	{
    	p = (*treenode)[i];
    	q = p->next;
   		while(q != p)
    	{
        	node * r = q;
        	q = q->next;
        	free(r);
    	}
    	free(p);
  	}
  	free(*treenode);
}
/**
 * Used in fitch 
 */
void freew(int nonodes, pointptr treenode)
{
  	int i, j;
  	node *p;

  	for (i = 0; i < (spp); i++) 
	{
    	free(treenode[i]->w);
  	}
  	for (i = spp; i < nonodes; i++) 
	{
    	p = treenode[i];
    	for (j = 1; j <= 3; j++) 
		{
      		free(p->w);
      		p = p->next;
    	}
  	}
}
/**
 * Used in fitch 
 */
void freed(int nonodes, pointptr treenode)
{
  	int i, j;
  	node *p;

  	for (i = 0; i < (spp); i++) 
	{
    	free(treenode[i]->d);
  	}
  	for (i = spp; i < nonodes; i++) 
	{
    	p = treenode[i];
    	for (j = 1; j <= 3; j++) 
		{
      		free(p->d);
      		p = p->next;
    	}
  	}
}
/**
 * Used in fitch and kitch
 */
void allocw(int nonodes, pointptr treenode)
{
	int i, j;
  	node *p;

  	for (i = 0; i < (spp); i++) 
	{
		treenode[i]->w = (vector)malloc(nonodes*sizeof(double));
		if ( treenode[i]->w == NULL )
		{
			syslog(LOG_ERR,"allocw: failed to allocate memory(1)\n");
			exit(-1);
		}
  	}
  	for (i = spp; i < nonodes; i++) 
	{
    	p = treenode[i];
    	for (j = 1; j <= 3; j++) 
		{
      		p->w = (vector)malloc(nonodes*sizeof(double));
			if ( p->w == NULL )
			{
				syslog(LOG_ERR,"allocw: failed to allocate memory(2)\n");
				exit(-1);
			}
      		p = p->next;
    	}
  	}
}
/**
 * Used in fitch and kitch
 */
void allocd(int nonodes, pointptr treenode)
{
	int i, j;
  	node *p;

  	for (i = 0; i < (spp); i++) 
	{
    	treenode[i]->d = (vector)malloc(nonodes*sizeof(double));
		if ( treenode[i]->d == NULL )
		{
			syslog(LOG_ERR,"allocd: failed to allocate memory(1)\n");
			exit(-1);
		}
  	}
  	for (i = spp; i < nonodes; i++) 
	{
    	p = treenode[i];
    	for (j = 1; j <= 3; j++) 
		{
      		p->d = (vector)malloc(nonodes*sizeof(double));
			if ( p->d == NULL )
			{
				syslog(LOG_ERR,"allocd: failed to allocate memory(2)\n");
				exit(-1);
			}
      		p = p->next;
    	}
  	}
}
/** 
 * Allocate spp tips and (nonodes - spp) forks, each containing three
 * nodes. Fill in treenode where 0..spp-1 are pointers to tip nodes, and
 * spp..nonodes-1 are pointers to one node in each fork. 
 * Used in fitch, kitsch, neighbor 
 */
void alloctree(pointptr *treenode, int nonodes)
{
	int i, j;
  	node *p, *q;

  	*treenode = (pointptr)malloc(nonodes*sizeof(node *));
	if ( *treenode == NULL )
	{
		syslog(LOG_ERR,"alloctree: failed to allocate memory(1)\n");
		exit(-1);
	}
  	for (i = 0; i < spp; i++)
	{
    	(*treenode)[i] = (node *)malloc(sizeof(node));
		if ( (*treenode)[i] == NULL )
		{
			syslog(LOG_ERR,"alloctree: failed to allocate memory(2)\n");
			exit(-1);
		}
	}
  	for (i = spp; i < nonodes; i++) 
	{
    	q = NULL;
		for (j = 1; j <= 3; j++) 
		{
		  	p = (node *)malloc(sizeof(node));
			if ( p == NULL )
			{
				syslog(LOG_ERR,"alloctree: failed to allocate memory(3)\n");
				exit(-1);
			}
		  	p->next = q;
		  	q = p;
		}
		p->next->next->next = p;
		(*treenode)[i] = p;
	}	
}
/** 
 * Initialize a tree. Used in fitch, kitsch, & neighbor 
 */
void setuptree(tree *a, int nonodes)
{
  	int i=0;
  	node *p;

  	for (i = 1; i <= nonodes; i++) 
	{
		a->nodep[i-1]->back = NULL;
		a->nodep[i-1]->tip = (i <= spp);
		a->nodep[i-1]->iter = true;
		a->nodep[i-1]->index = i;
		a->nodep[i-1]->t = 0.0;
		a->nodep[i-1]->sametime = false;
		a->nodep[i-1]->v = 0.0;
		if (i > spp) 
		{
			p = a->nodep[i-1]->next;
			while (p != a->nodep[i-1]) 
			{
				p->back = NULL;
				p->tip = false;
				p->iter = true;
				p->index = i;
				p->t = 0.0;
				p->sametime = false;
				p = p->next;
			}
		}
  	}
  	a->likelihood = -1.0;
  	a->start = a->nodep[0];
  	a->root = NULL;
} 
/**  
 * Read in distance matrix. Used in fitch & neighbor 
 */
void inputdata(boolean replicates, boolean lower,
    boolean upper, vector *x, intvector *reps)
{
	int i=0, j=0, columns=0;
	boolean skipit=false, skipother=false;

  	if (replicates)
    	columns = 4;
  	else
    	columns = 6;
	for (i = 0; i < spp; i++) 
	{
    	x[i][i] = 0.0;
    	scan_eoln(infile);
    	initname(i);
    	for (j = 0; j < spp; j++) 
		{
      		skipit = ((lower && j + 1 >= i + 1) || (upper && j + 1 <= i + 1));
      		skipother = ((lower && i + 1 >= j + 1) || (upper && i + 1 <= j + 1));
      		if (!skipit) 
			{
        		if (eoln(infile))
          			scan_eoln(infile);
        		if (fscanf(infile, "%lf", &x[i][j]) != 1)  
				{
          			syslog(LOG_ERR,"The infile is of the wrong type\n");
          			exit(-1);
        		}
        		if (replicates) 
				{
          			if (eoln(infile))
            			scan_eoln(infile);
          			if (fscanf(infile, "%d", &reps[i][j]) != 1) 
					{
            			syslog(LOG_ERR,"The infile is of the wrong type\n");
            			exit(-1);
          			}
        		} 
				else
          			reps[i][j] = 1;
      		}
      		if (!skipit && skipother) 
			{
          		x[j][i] = x[i][j];
          		reps[j][i] = reps[i][j];
      		}
      		if ((i == j) && (fabs(x[i][j]) > 0.000000001)) 
			{
       			syslog(LOG_ERR,"diagonal element of row %d of distance matrix "
        		"is not zero.\n"
        		"       Is it a distance matrix?\n",i+1);
        		exit(-1);        
      		}
      		if ((j < i) && (fabs(x[i][j]-x[j][i]) > 0.000000001)) 
			{
        		syslog(LOG_ERR,"distance matrix is not symmetric:\n"
        		"       (%d,%d) element and (%d,%d) element are unequal.\n"
				"       They are %10.6f and %10.6f, respectively.\n"
				"       Is it a distance matrix?\n",i+1, j+1, j+1, i+1,x[i][j], 
					x[j][i]);
        		exit(-1);
      		}
    	}
  	}
  	scan_eoln(infile);
} 

